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Abstract. By applying effective medium— style calculations to random spring 
networks, we demonstrate that internal stresses fundamentally alter the nature 
of the rigidity transition in disordered materials, changing it from continuous 
to first-order and increasing the mean coordination number z at which rigidity 
first occurs. Furthermore, we predict the existence of a novel stability regime at 
low z when the distribution of stresses is asymmetric. Means of verifying these 
predictions are suggested. 



PACS numbers: 46.32,+x, 61.43.-j, 62.20.Fe 

1. Introduction 

Predicting under what conditions a given material will support an applied load without 
undergoing plastic deformation is clearly of great importance to materials science 
and industry alike. When the material in question is ordered, the periodicity of 
the microstructure allows the elastic moduli to be derived from the properties of 
the constituent particles in a manner that is, at least in principle, straightforward. 
Inhomogenous materials pose significant new problems p^. Model disordered systems 
have revealed that the elastic moduli typically vanish continuously as the volume 
fraction (or some related parameter) drops below a critical value |2J. This rigidity 
transition appears to bear some of the hallmarks of continuous phase transitions in 
thermal systems [5] , such as a diverging fluctuation correlation length and some degree 
of universality in scaling behaviour near the transition 4 . Experiments on materials 
such as chalcogenide glasses are broadly consistent with this picture Ej (but see 
later) . 

In all of the model systems considered thus far, it has been assumed that the 
material is initially either everywhere unstressed, or everywhere under tension 7 . 
An overlooked possibility is the existence of internal stresses, i.e. stresses that exist 
on the microscopic scale that average to zero on macroscopic lengths. This omission 
would be valid if the material was formed in a way that allowed the total elastic 
energy stored in interparticle bonds to fully relax to its global minimum. However, 
there has been much recent interest in materials for which this is unlikely to be true. 
Glasses [8], particulate constructs [§], and soft condensed matter systems such as 
dense foams and emulsions ^U] are formed by the kinetic arrest of the constituent 
particles under suitable driving and boundary conditions, producing configurations 
which have a quenched-in distribution of interparticle forces. Indeed, simulations of 
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ground-state atomic glasses have revealed internal stresses far greater in magnitude 
than the macroscopic value and the long-time relaxation of a range of 'jammed' 
soft materials has been generically attributed to stressed local regions 

That internal stresses will alter the mechanical response is clear. Elastic filaments 
provide a familiar and canonical example: when placed under tension they may be 
plucked and will return to their original state, in contrast to compressed bands which 
readily buckle Such considerations becomes no less important on microscopic 

length scales, where internal forces can significantly alter the mechanical stability of 
the ensemble, even in the absence of a macroscopic prestress m] . States of self-stress 
are already known to alter the classic Maxwell counting rules, increasing the number 
of displacement modes that do not violate the constraints imposed by interactions |15| . 
However, it is not yet known how internal stresses alter the nature of the transition. 

The purpose of this paper is to determine the role of internal stresses on 
the rigidity transition to and from mechanically stable states for a simple class of 
disordered solid, namely random Hookean spring networks at zero temperature. This is 
achieved using a novel, analytically-tractable approximation scheme that qualitatively 
agrees with known results for unstressed systems. We find that any internal stress 
changes the transition from continuous to first order, i. e. the elastic moduli are finite 
at the onset of rigidity. Furthermore, the mean number of contacts per particle z at 
which the system becomes rigid increases with respect to its unstressed value, with 
a power law dependency on the magnitude of internal stresses. We also predict the 
existence of a new stable regime at low z for asymmetric stress distributions. Means 
to experimentally verify these predictions are suggested. 

2. Assumptions of the model 

Given that the approach adopted in this article involves several novel approximations, 
it is perhaps sensible to first discuss their strengths and limitations with respect to 
more established approaches. These will be discussed in this section. For future 
reference, the combined approximation scheme used here shall be referred to as the 
mean mode approximation or MMA. 

The closest scheme to MMA is known as the effective medium approximation or 
EMA. Here, the inhomogeneous network is mapped onto an analogous, homogeneous 
substrate for which the Green's response to a point force is known. For bond-diluted 
lattice models, the choice of analogous system is clear: a complete lattice with every 
bond occupied, but where every spring has an effective spring constant k eS . The goal 
is then to find k eS . A key problem in extending this method to prestressed systems 
is that the Green's function is usually not known, even for a homogeneous system [7] . 
The MMA approach, on the other hand, instead assumes that the displacement field 
around a perturbed bead takes a particular, averaged form, as detailed below. This 
can be viewed as a more drastic approximation than EMA, and indeed the reduction 
in the local degrees of freedom inevitably results in the rigidity transition occuring 
at a lower coordination number z as EMA. However, it has the advantage of not 
requiring an analogous system (with a known Green's function) to be identified, and 
can therefore be applied to unstressed and internally stressed networks alike. For 
unstressed networks, for which results from EMA are known, the MMA approach 
adopted here gives qualitatively the same behaviour as EMA, including a transition 
at finite coordination number z and an exponent of unity for the elastic moduli as 
one traverses the transition. It is therefore not unreasonable to suppose that it also 



First order rigidity transition in internally stressed networks 



3 



qualitatively predicts the correct behaviour when internal stresses are incorporated. 

Much of our understanding of rigidity percolation has come from simulations, so 
it is natural to ask why a numerical scheme has not been adopted. The answer is 
essentially one of simplicity. For unstressed networks, it is straightforward to simply 
dilute bonds from an ordered lattice, taking care to slightly displace the lattice nodes 
to avoid colinearity and coplanarity of bonds. Force balance is obeyed everywhere by 
the simple fact that all forces are zero. This is not true in the presence of stresses: 
removing a bond from a stressed network will typically violate local force balance, 
requiring relaxation of the network to a new stable configuration for every contact 
broken. Thus networks constructed by a physically realistic procedure that do not 
relax all of their internal stresses are fundamentally a product of their history, and 
CPU-intensive algorithms such as molecular dynamics will be required to generate 
the proper starting condition. Such methods are beyond the scope of this article. The 
aim of this paper is partly to inspire simulation work to verify the predictions made, 
but also to pave the way for a full, theoretical understanding of the rigidity transition 
in stressed systems, for which the MMA provides a first, quite possibly qualitatively 
correct, first step. 

Finally, we comment on the use of central force networks. The history of rigidity 
percolation has been somewhat confused by the study of central force networks, 
although the situation for unstressed networks is now clear (see e.g. pQ for a full 
range of results). An obvious objection is then, why risk jeopardising this work with 
a potentially troublesome system? The answer is twofold. Firstly, simplicity: central 
force systems have fewer parameters to consider and therefore give a clearer picture 
of the effects of internal stresses. Secondly, however, and much more importantly, is 
that many of the systems mentioned in the introduction (for which internal stresses 
are likely to exist), such as wet foams, emulsions, frictionless granular media and 
Lennard- Jones systems |SI EI , are central force. Studying central force systems 
is therefore of immediate and significant interest if one hopes to understand these 
materials. It also serves as a first step towards understanding other materials which 
are not central force, such as frictional granular media and atomic systems, which 
could be tackled by constructing an extending MMA scheme with transverse forces. 

3. Definition of the model 

Our model system consists of a collection of particles a interacting via some known 
finite-ranged interaction potential to produce a static body in which force balance 
is everywhere obeyed. (In the language of networks, the a are nodes interconnected 
by force-transferring bonds.) The position of each particle is specified by the d- 
dimensional Cartesian vector x Q , and the displacement to a connected particle j3 
is x' 3 — x Q = r Q/3 n a/3 in terms of the unit vector and the centre-to-centre 
distance r a/3 > 0. The force on (3 due to a is given by the central law f(r) as 
fa/3 _ j( r Q / 3 ) n Q / 3 ; so compressive contacts correspond to positive /. For simplicity, 
we assume all interactions are Hookean with identical spring constants k > and 
natural lengths ro, i.e. f(r) = —k{r — Tq) = —kr$e in terms of the dimensionless 
extension e = (r — ro)/Vo. 

The macroscopic stress in the initial configuration depends on the joint probability 
distribution of h a ^ and r Q/3 . For simplicity, we assume here that the bond vectors 
n"^ are uniformly distributed on the unit (d — 1) -dimensional hypersphere, and 
independent of the r Q/3 . This ensures the macroscopic stress field is isotropic. 
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Since our ultimate goal is to elucidate the role of internal stresses, we look for the 
simplest distribution P(e a/3 ) that gives zero global stress and thus choose P{e) — 
pS(e — £o) + (1 — p)5(e — £1) with p 6 (0, 1). The macroscopic stress tensor is then 
diagonal with pressure ~ kr^~ d [peQ + (1 — p)e{\ for <C 1, and requiring that 
this vanishes gives e% « — Eqp/(1 — p) in terms of the two free parameters £o and p. 
Note that here and below we consider \ei\ <C 1, i.e. a small perturbation around the 
unstressed limit Si = 0, as this facilitates some of the subsequent analysis, although no 
crucial modification for larger deformations are expected. Physically, this corresponds 
to systems of slightly deformed particles, or to a network that has relaxed close to, 
but not quite reaching, its global energy minimum in which all nodes are separated 
by their natural spring lengths. 

The mechanical stability of a configuration is determined by applying an 
infinitesimal external force 5f cxt to a randomly selected particle a, and allowing all 
particles (3 (including a) to relax into a nearby configuration + dxP . If this nearby 
configuration obeys force balance, and the work done by the external agent 5W is 
positive, the system is deemed stable; negative work is assumed to signify mechanical 
instability due to the inability of the interparticle forces to support the load. Note that 
this force probe need not come from some spontaneous fluctuation in the steady state, 
whose origin would be obscure for the athermal systems under consideration here, 
but may be due to mechanical noise from an external source, or the final fluctuation 
in contact forces as the system is quenched into its final configuration by whatever 
preparation procedure is employed. 

In analogy with statistical mechanics, the macroscopic response is expected 
to depend on the microscopic configuration only through a small set of suitably 
selected parameters. We henceforth ensemble average over all local degrees of freedom 
consistent with a given bond extension distribution P(e) and the mean coordination 
number (number of bonds per particle) z. The requirement of force balance in the 
perturbed configuration is 

where the sum is over all [3 interacting with a. Here, the change in interaction force 
on (3 due to a has been written as = —kA^(5x^ — Sx"), where for the central 

forces under consideration here, 

At = nfnf + e a0 (f« - nfnf) , (2) 

with Sij the Kronecker delta. This projects out the component of the relative 
displacement parallel to the bond, which alters the magnitude of the force by an 
amount proportial to — k, and the transverse component, which (in this linearised 
scheme) describes the rotation of the original force — kr^e 01 ^ . 

The displacement field <5x depends on the entire initial configuration and is too 
complex to treat exactly. Instead, observe that (Sx a ) = A<5f oxt for the isotropic 
systems under consideration here, where the value of the compliance A is determined 
later. A first approximation is then to replace the global dependency of <5x Q in Q 
by the local form <5x" = XSf cxt . No such average form is immediately apparent for 
£x^, but closure is possible by assuming that the change in contact force 5f al3 can 
be viewed as an external force on /3, i.e. Sx 13 = A<5f Q/3 with the same A as before. 
Inserting these two approximations into and (J2J, and eliminating the <5x's gives 
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sff = s;y,\r/'. 

sf = [i + ixky^nfnf 

+ [l + (e^\k)-Y 1 (s t] -nfnf) (3) 

The second rank tensor S^J 3 gives the propagation of force from a to a connected 
particle /3, and is expressed here in terms of projection operators parallel and 
perpendicular to the contact vector h a ^ . Note that the unphysical singularity at 
sXk = — 1 is avoided by the stability equation given below. 

Inserting Q into (JIJ and averaging over P(s) and n a/3 gives the following equation 
for A, 

The work done by the external agent is 8W = ^A|5f ext | 2 , so stability corresponds 
to positive, real roots of A in Although A is in principle a measurable quantity, 
a more convenient order parameter for quantifying the order of the transition is the 
bulk elastic modulus K, which is related to A via K ~ rQ~ d \~~ 1 with a dimensionless 
prefactor that will not concern us here. 



4. Results 



To test the validity of the approximations, we first consider unstressed systems 
£q = £i = 0, for which the transition is known to be continuous. In this case, (@J 
reduces to a linear equation with the single solution Afc = (z/d— 1) , admitting stable 
systems only when z > d. The corresponding modulus is K ~ kr^~ 2 {z/d — 1), which 
can be written as a power law K ~ (z — z c y to explicitly demonstrate the continuous 
transition at z = z c = d with an exponent / = 1. Effective medium theory predicts 
the same exponent but at the higher transition point z° ma = 2d 16 1, as predicted 
by Maxwell counting. However, EMA requires knowledge of the Green's function for 
an equivalent homogeneous system. The advantage of our MMA approach is that it 
predicts a fmite-z transition without requiring a known Green's response, and thus 
can be applied to a broader range of materials, at least qualitatively. In particular, 
the generalization to internal stresses is straightforward, as we now demonstrate. 

The simplest state of internal stress to consider is a symmetric distribution 
So = — £i, i.e. p = |. In this case, equation only permits real positive roots 
in one region of parameter space corresponding to z > z c (eq) > d as plotted in Fig.^ 
The transition points lie along the curve 

4(2 ~ d)3 (5) 



u (3d) 3 (^-l)' 

which obeys the expected symmetry under £o «-» — £o- The value of K at the transition 
is most easily expressed in terms of z — d, 

is ji.-d\-\ 2fcrp d (z — d) 

-Ktrans ~ r Q A = — (bj 

6d 
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Figure 1. Combinations of coordination number z and dimensionless extension 
£0 that generate stable systems, here shown for the symmetric case V = \ ( so 
ei = — so) and dimension d = 3. The solid disc is the continuous transition; all 
other transitions are first order. (Inset) K versus z for eo = (solid straight line) 
and so = 0-01 (dashed), demonstrating the first-order transition and multiplicity 
of solutions in the latter case. 



By inspection of and ©, it is clear that the transition is first order everywhere 
except at eo = 0, i.e. internal stresses both move the transition to a higher z c , and 
change its nature to first order. 

This result broadly holds for asymmetric distributions with p ^ |, i.e. e ^ 
although the transitions for eo > and eo < are no longer symmetrical. Expanding 
about z — d to 0[(z — d) 2 ], the transitions are at 



= / 4(l-p)(z-d)3 {l-2p)(z-df 
V p{MY{d-l) 3pd 2 {d-l) {> 

The leading order term ~ (z — <i) 3 / 2 admits two roots, corresponding to the transition 
lines for eo > and eo < 0. The next-to-leading order term always has the same sign 
(that of 1 — 2p) for both roots, thus breaking the symmetry about eo = 0. K at the 
transition is also asymmetric, but remains always first -order for Co 7^ 0, 



2krl- d {z-d) f l-2p 3(z-d) \ 

K t ^ s x|l + ^ ^—^—^ I ( 8 ) 

plus higher order terms, where the sign of the square root is chosen to match that in 

o 

However, the modulation of the transition curves for z > d is not the only effect 
of asymmetric stresses; a distinct stable regime with z < d also emerges, as seen in 
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Figure 2. As Fig. but for an asymmetric distribution of internal stresses in 
which a fraction p = 0.1 bonds have dimcnsionlcss extension so an d 1 — p = 0.9 
have ei = — eop/(l — p)- The thick line and solid disc correspond to continuous 
transitions; all other transitions are discrete. 



Fig. [3 The £o = boundary of this region stretches from z = 1 to z = z* with 

4<2p(l - p)~ 



l = (d-l) 



1 



(1 - 2p)2 



(9) 



and extends strictly in the direction of sign opposite to 1 — 2p (i.e. eq < for p < i), 
narrowing as Eq increases in magnitude. In words, this regime corresponds to systems 
with a small proportion of highly compressed contacts in a sea of weakly tensile bonds 
(i.e. small p and e < 0). We speculate that the z > d and z < d stable regions 
may be analogous to the glass and gel states, respectively, in colloids with short-range 
attractions |17j , although the transition here is a purely percolation phenomenom and 
has no entropic component. That this new regime extends down to the unphysical 
value z = 1 is most likely a further consequence of the approximations involved: just 
as the unstressed transition lies at z = 2d in real systems, so would the lowest allowed 
value be z = 2, as realised in e.g. long strings of beads under tension. 

5. Discussion and summary 

Experimentally verifying these predictions should, in principle, be straightforward. 
According to this theory, the magnitude of K at the transition, and the increase in 
the transition value z c , should both increase monotonically with the magnitude of 
the internal stresses. Indeed, this effect may have already been observed in Si^Sei-a; 
glasses [5], although it is currently interpreted as evidence for an intermediate rigid- 
but-stressless regime resulting from self-organization of the contact topology during 
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cooling JSj. Differentiating between these two descriptions could be achieved by 
varying the rate of quench [6]: faster quenches will generate a broader distribution 
of internal stresses, and hence (according to the theory presented here) a higher z c 
and greater jump in elastic properties at the transition. Slow quenches would give 
small internal stresses and could even appear as continuous transitions, to within 
experimental error. It is less clear how the degree of asymmetry p ^ \ wm depend on 
sample preparation, and hence it is difficult to predict when the low-z stable regime 
will occur. Molecular dynamics simulations, where some quantity analogous to p could 
be directly measured, may help to resolve this issue. 

Unless there exist residual interactions capable of maintaining solidity, systems 
placed in an unstable region will plastically rearrange according to the dynamical 
properties of the particles and the surrounding medium. Clearly such dynamics 
cannot be described by the static theory presented here. Nonetheless, some qualitative 
observations based on our results can be made. Firstly, unless the system completely 
relaxes its internal stresses during the dynamical phase, it will not have a coordination 
number at or near to the usually quoted percolation transition (z = 2c? here), since such 
systems are simply not stable when there are internal stresses. This is a qualitative 
prediction that could be verified by molecular dynamics simulations, for instance. 
Secondly, assuming that the rearranging system eventually comes to halt on the 
boundary between stable and unstable regions, as recently proposed for the clustering 
of weakly-attractive colloids , then the elastic moduli will only become arbitrarily 
small if the magnitude of internal stresses are similarly small. Again, this is a clear 
prediction of a qualitative difference between stressed and unstressed networks. 

A secondary aspect of this article has been the introduction of a new 
approximation method, the MMA, which predicts a range of non-trivial behaviour 
despite the simplicity of its assumptions. Indeed, it qualitatively reproduces the 
results of effective medium theory for unstressed systems with much less algebra. 
It is therefore sensible to ask when the approximation is expected to work, and when 
it might fail. The MMA proceeds by closing the force balance equations under the 
assumption of a parameterised form of local response. While this will always fail 
quantitatively, it should only qualitatively fail when the actual response is very different 
to the assumed form. In this instance something similar to the actual mode (if known) 
could be employed instead. Another potential problem is that the displacement of 
particles connected to the perturbed one is assumed to depend only on the change 
in contact force with the perturbed particle itself. Given that the mean force must 
decay monotonically with distance (to obey force balance), the closest bonds will be 
perturbed the greatest and so this seems reasonable; however, exotic modes in which 
the force does not decay in every direction may cause this assumption to fail. Even in 
such instances, it is hoped the MMA could be extended to faithfully mimic the actual 
response. 

In summary, we have argued that internal stresses qualitatively alter the nature 
of the rigidity transition to configurations with non-vanishing elastic modulii, making 
it first order and also moving the threshold coordination number z c to a higher 
value. A distinct stability regime with low z was also predicted to arise when the 
internal stresses are asymmetrically distributed. Although only central forces have 
been considered here, these basic findings are expected to extend to systems with 
bending forces and other microscopic interactions. It is hoped that numerical and 
experimental verification of these claims will be forthcoming. 
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